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Abstract 

We study a one-dimensional lattice flocking model incorporating all three of the 
flocking criteria proposed by Reynolds [Computer Graphics 21 4 (1987)]: alignment, 
centring and separation. The model generalises that introduced by O. J. O' Loan 
and M. R. Evans J. Phys. A. 32 L99 (1999). We motivate the dynamical rules by 
microscopic sampling considerations. The model exhibits various flocking regimes: the 
alternating flock, the homogeneous flock and dipole structures. We investigate these 
regimes numerically and within a continuum mean-field theory. 

Pacs: 05.70.Fh 64.60.Cn 87.10e 

1 Introduction 

Flocking refers to a family of behaviours regularly seen in nature, including schooling of 
fish[T] and nocking of birds and moths[2 4 . Whilst the features that are truly essential to 
flocks remain open to speculation [3J EJ |H| it is clear that the variety of biological[7l [TJ |2] 
and engineered [H] systems in which flocking behaviour does, or could, prove useful is vast. 
Statistical mechanics has recently played an important role in uncovering and determining 
the source of many flocking phenomena jH|. 

One of the pioneers in flocking research was Reynolds 0. Within his work he determined 
three behaviours characteristic of flocking: the desire to match the velocity of flockmates 
(alignment), the desire to stay close to flockmates (centring), and the desire to avoid colli- 
sions (separation). Although Reynolds' original research was motivated by applications in 
computer graphics, his definitions have been adopted widely and implemented by researchers 
from a variety of disciplines. Reynolds coined the term boids to describe the generic self- 
propelled particles obeying these rules, and we shall continue that convention in this paper. 

Within the physical science literature research has focused primarily on the steady state 
behaviour of flocking models and the non-equilibrium elements which cause solutions to differ 
from equilibrium systems. Many models involving only an alignment interaction have been 
shown to demonstrate a symmetry broken velocity statejTH]. Results have also been found 
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which are strongly non-equilibrium in one and two dimensions, where equilibrium systems 
do not undergo comparable symmetry breaking phase transitions^] • Important work on 
alignment has included that of Vicsek et al^U] and Toner and Tu^T] who shed light on 
the non-equilibrium processes involved. A wide variety of models also exist in two or three 
dimensions incorporating one or several of Reynolds' rules besides alignment |13l 114*]. The 
lattice model of which includes separation and alignment, shows a density-dependent 
symmetry broken state, for example. Similar results have also been derived using a network 
approach to flocking, whereby boids are nodes within a dynamic graph [To]. 

There are fewer obvious physical realisations of one- dimensional flocking, although anal- 
ogy may be made with ring shaped aquariums, and with pedestrian dynamics in corridors 
[IB] . However, research into one- dimensional systems certainly is interesting from a funda- 
mental viewpoint: in particular, nonequilibrium phase transitions may be exhibited|17j. In 
the one dimension alignment model of Czirok et al^B] and the one dimensional lattice model 
of O'Loan and Evans [T^] novel states not anticipated by higher dimensional treatments were 
found. In [TH] a symmetry broken steady state was reported, whilst [19 a found a steady state 
which alternated its direction of motion on an unusually short time scale O(lniV) where N 
is the number of boids in the system. In the one dimensional model of Levine et al[2I] the 
effects of centring and separation were considered and new behaviours were also found. 

Our aim in this work is to introduce and study a one-dimensional lattice model involving 
all three of Reynolds' flocking criteria. We do this by generalising the model of jTU]. We 
demonstrate a variety of new flocking regimes: further to the alternating flock state found 
in pH] we show that on small systems a homogeneous flock may occur; also, the inclusion 
of centring may produce various dipole structures in which the velocities of the boids point 
towards the centre of the structure. Our generalisation also addresses some deficiencies in 
the dynamics of ^H] which implied full knowledge of neighbouring boids — in the present 
work the dynamics is motivated from the sampling of a finite number of neighbouring boids 
and yields a smoother form [Eq. © below] for the expected update velocity of a boid as a 
function of the average neighbourhood velocity. 

The paper is organised as follows. In section 2 we define and motivate the model to 
be studied. In section 3 we present numerical simulations of the model and examine the 
different regimes that may be observed. In section 4 we employ mean-field theory to provide 
further evidence for these regimes. We conclude with a discussion in section 5. 

2 Boid Dynamics 
2.1 Model synopsis 

We first review the model of JH] which consists of iV boids inhabiting a lattice of L sites 
with unit lattice spacing. The global density, (p), is defined as (p) = N/L. Each boid /x is 
defined by a position — 1 . . . L and velocity = ±1. An update of the system consists 
of the following steps: 

1. A boid is chosen at random. 
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2. A preferred direction U(x^,t) = ±1 or is determined from the velocity and spatial 
distribution of neighbouring boids. 

3. Vfj, is updated to ±1 with probability W±i 

W ±1 = l/2±{l/2-r } )U(x li ,t) (1) 
where the constant r\ is a parameter of the model. 

4. The boid moves: x^ is updated to x^ + t> M . 

We take the neighbourhood of site i to be the sites i — 1, i and i + 1. In the model of ^H] 
U is given by the local majority velocity which we now define. Let the number of boids on 
a site with velocity v be n v (x,t), the site density p(x,t) = ^2 v n v (x,t) an d site momentum 
<f>(x, t) = vn v (x, t). The majority velocity of neighbouring boids is equivalent to the sign 
of the neighbourhood (3 site) averaged momentum <p(x, t) given by 

1 1 

(f>{x,t) = - ^ <f>(x + ia,t) . (2) 
»=-i 

In what follows such a neighbourhood average of a quantity X is represented by the hat 
symbol X. The preferred direction U(x,t) is then determined by 

{1, if0O,t)>O 
-1, if?(x,t)<0 (3) 
0, if?(x,t) = 

Thus in the model of ^H] the preferred direction is given deterministically by the local 
momentum and the model thus includes alignment. Following a determination of U, rj is 
the probability that the velocity is not updated to this value (and instead to —U). This 
probability is independent of the dynamical variables. 

In this paper we generalise the model of [1.9.] at step (ii) above to include the effects of 
centring and separation as well as alignment — this will be discussed in section I2~4l We also 
generalise to consider the case where U is itself a stochastic function of local variables: a 
boid may determine its preferred direction from the local variables as ±1 or with some 
probabilities (as will be motivated below). However, since U appears linearly in (JIJ), averaging 
over its possible values amounts to replacing U by its expectation value U in ([TJ: 

W ±l = l -[l±G{x^t)] (4) 

where 

G(x,t) = (I-277) U. (5) 

Thus the important quantity is U{x, t), the expectation value of the preferred direction which 
should be a function of the local density and velocity, for example. The expectation value 
of the updated velocity is given by G(x,t). 
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In the next subsection we propose an explicit form for U and hence G(x,t). First it is 
useful to review the form of (|5Jl. The quantity sgin(G) is the velocity to which a boid is 
updated in the absence of noise, whereas \G\ is the certainty with which this outcome is 
attained. The uncertainty in the outcome comes from two sources, the first is due to random 
errors, controlled by the parameter rj, where the boid moves against its preferred direction U. 
The second is from the stochastic nature of U itself, which as we argue in section 12.31 comes 
from the uncertainty with which a boid perceives the local flock properties and therefore 
determines its preferred direction. 



2.2 Alignment 

In many models fHl El E] it has been shown that alignment alone is sufficient to pro- 
duce complicated behaviour recognisable as flocking. Alignment is effected in our dynamics 
through G. In this work we consider G to be of the form 

G = (l-2,) tanh( ^ t)) (6) 
v /; tanh(/3) v 1 

where 

p(x,t) 

is the neighbourhood average boid velocity. The parameter (3 > controls how nonlinear 
G(x,t) is as the function of V(x, t). In the limit where (3 — > oo the majority rule case © 
is in effect. In the limit (3 — > we obtain a linear function. The form (JBJ) addresses several 
deficiencies in (j3J). Firstly © is analytic, whereas J3J) is non-analytic. Secondly with (0) 
boids become sensitive to increasing majority size, as might be expected in physical systems. 
In the following section we will see how this second intuitive feature can arise out of errors 
in local flock perception. 



2.3 Justification of form of G via sampling argument 

We now turn to a justification of the form © from microscopic considerations. We present 
an argument amounting to a simple algorithm for the determination of G based on sampling 
of neighbours. In the model of fTS] boids determine U with perfect knowledge of all neigh- 
bouring boids (and perfect ignorance of other boids); hence a local majority is determined 
with certainty. Consider instead the situation where a boid can only make M observations 
of its neighbours, each observation is of one boid selected randomly from neighbouring boids 
with replacement. (Also note that from the definition of neighbourhood a boid acts as a 
neighbour to itself.) The Majority rule applied to the sample of M boids then requires a 
fixed observation time O(M) regardless of local density. This is in contrast to the strict 
majority rule algorithm which requires information about all neighbours 0, there being no 
restriction on the potential number of neighbours. 
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Consider a sample of M neighbouring boids from the group of 3p(x, t) neighbours, con- 
taining 3ni(x,£) rightward and 3n_i(x,t) leftward moving boids. With binomial probability 
P(k), k will be selected from the rightward and M — k from the leftward moving groups: 

™ -CO (!)'(¥)" 

The preferred direction is determined from the sample as Um = sign(2/c — M). The expected 
updated velocity Gm is, taking M is odd, 

M-1 
2 

G M = (1 - 2 V )U M = (1 - 2tj) ^ [P(M - k) - P{k)\ . (9) 

k=0 

Using the substitution 

^ 1+jVQM) 
%(a;,t) = p(x,t) , 

Gm can be expanded as an odd polynomial in V: 

M-1 M-1 k 

k=0 i=0 j=0 \ / \ / \ J ' / 

A comparison of the two forms of G ()10|) and © shows a qualitative correspondence 
between increasing M and f3. For example, we have from (|1(J|) 

Gi(x,t) = (l- 277)y(x,t) 

G s (x,t) = (l-2^)^|^(3-y(x,t) 2 ) 

GfiCM) = (1-2^)^^(15- W(x,t) 2 + 3V(x,t) 4 ), 

o 

whereas expanding (jHJ) in powers of V" yields 

G(z, t) = (1 - 2^)^^(1 - ^ + l^*' ^ + *)")) • 

In the two limiting cases Gm^oo = hnig^oo G(x, t) and Gm=i = linig^o G(x, t) the corre- 
spondence is exact. 

For intermediate values of M it is possible to define an approximate correspondence 
through a function (3{M) using one of several schemes. Anticipating our mean field treatment 
below (section 0J) we choose to match the gradient of G with respect to v at the origin: 

M-1 

k=0 V ' \ 2 / 

This yields, for example, /3(M = 3 or 4) = 1.5 and (3{M = 25 or 26) « 4. 

To summarise we have argued that the form (JHJ), parametrised by /3, for the expected 
updated velocity qualitatively corresponds to (but has a simpler form than) the result of 
taking a sample of M neighbours. The correspondence (3{M) is quantified in 
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2.4 Centring and Separation 



Centring, the tendency to move towards the local centroid u(x,t), defined as 



u(x, t) 



can be introduced through G in an analogous way to alignment. If a relative importance 
of k is assigned to the effect of centring over alignment in any observation then (JHJ) can be 
generalised to 

G(x,t) = {1 _ 2v) ^mi-«W,t) + Ku,(z,t))) _ (12) 

The desire for separation, typically through hard core repulsion, is very restrictive in one 
dimension. If instead we consider separation as some cost associated with moving through 
a dense region then it amounts to a tendency to move away from the centroid (opposite to 
centring). In order to include these effects together we introduce a capacity, C, which is 
the capacity at which the relative strengths of these two effects neutralise. When the local 
density t) is greater than C the desire for separation exceeds the desire to centre, and 
vice versa. Thus, our most general form for G is chosen as 



G(x,t) = ^ , 2 ^ tanh<i/3 
tanh p 



[l-K)V{x,t) + K C J^M x,t) 



(13) 



Note that is the limit C -> oo of (IT3J) . 

At this stage it is useful to summarise the parameters of the model with our choice of G: 
the capacity C is the density where the separation tendency surpasses the centring tendency; 
k determines the relative strengths of spatial density (centring and separation) and alignment 
effects. Thus k and C determine the relative importance of Reynolds' three flocking criteria. 
The parameters rj and j3 may both be considered as sources of noise acting through G: (3 
introduces sensitivity to different field strengths while rj introduces uncorrelated errors. With 
reference to section 1^73*1 it is useful to consider the noise from f3 as arising out of stochastically 
sampling the local flock to determine the boid's preferred direction whereas rj encompasses 
additional error sources once the preferred direction has been selected. 



3 Computer Simulations 

3.1 Observed Regimes 

Numerical simulations were undertaken for a variety of noise levels (/3 and rj), and flocking 
parameters (k and C). Run time constraints limited our investigations to systems of size 
upto 1024 sites and times upto 10 7 timesteps, where one timestep is sufficient time for each 
boid to be updated once on average (N updates). We found four characteristic and robust 
behaviours which we now discuss, alongside a myriad of intermediate behaviours. 
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Figure 1: Figures show the development of systems in space time plots from random initial 
density and velocity distributions on systems with L = 128 sites. The greyscale measures 
density p(x, t) on a logarithmic scale. In this and following figures parameter sets are indi- 
cated [thus] (a) [77 = 0.2,(3 — 1, k — 0.5, (p) — 1, C — > 00] A disordered system showing 
no sustained correlations in density or velocity, (b) [77 = 0.02, f3 — 8, k — 0.5, (p) = 1, 
C — > 00] An alternating flock exhibiting a non-zero velocity between regular reversals. 
(c)[?7 = 0.125, P — > 00, k — 0, (p) = 8] A homogeneous flock having homogeneous den- 
sity and fixed global velocity, emerging from transient alternations, (d) [77 = 0.02, f3 = 8, 
k = 0.75, (p) — 1, C — > 00] A dipole system consisting of several dipoles separated by low 
density homogeneous domains undergoes a slow coarsening process towards a single large 
dipole. (For the definition of a dipole see text.) 
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The first characteristic behaviour (figure (T^,) is the disordered state. This state has 
homogeneous density and a global velocity of zero, global velocity being the mean velocity 
of a boid in the system 

M = ®& . (14) 

This state persists at noisy parameter sets (high 77, low (3) especially where Reynolds' effects 
exist at equivalent strengths or where separation dominates (k > 0.5 and C < (p)). 

The alternating flock is the second regime (figured). This system undergoes a repeating 
pattern of cohesive traversals interspersed by rapid reversals, similar to that of [Hi] . 

The homogeneous flock as shown in figure ^ is a flocking state which complements the 
alternating flock. The homogeneous flock has a homogeneous density and fixed non-zero 
global velocity, (v). This regime is found to be unstable where noise is high or where global 
density is low, and is resilient only to small levels of centring or separation. The homogeneous 
flock is often established following an alternating transient as shown in the lower part of figure 
[^:. The homogeneous flock can also be observed in the model of but at larger densities, 
or smaller system sizes than those considered there 

The final regime observed is a dipole state (figure ^i), which exists in systems where 
centring is the dominating Reynolds' effect. A dipole is a dense, localised structure. As we 
shall discuss in detail (in section I3.3|) a dipole has a well defined centre with G pointing 
inwards trapping boids, and a sharply decaying density profile at its edges. A system of 
one or many dipoles can exist on the array, either as direct neighbours or separated by 
depopulated domains. As can be seen in (figure HJl) a slow coarsening process is observed in 
a system of many dipoles, whereby small dipoles are eliminated and the boids are absorbed 
by larger dipoles. The dipole states exist for large values of k and C i.e. relatively weak 
alignment and separation. 

These four regimes all have readily characterised behaviours and occupy characteristic 
regions of parameter space as shown in figure |2] which describes the case C — > 00 (absence 
of separation). Figure 121 may be thought of (and we shall refer to it) as a phase diagram. 
However, the different regimes are not truly "phases" as the boundaries between them depend 
on system size and, as we shall see below, the homogeneous flock regime actually disappears 
in the limit of large system size. In the following subsection we shall discuss in more detail 
the flocking and dipole regimes. 



3.2 Flocking regimes 

The alternating flock was originally studied in ^5] where only alignment was present. Here 
we find that it persists in systems with all three Reynolds' effects present provided alignment 
dominates (k < 0.5) and (3 remains finite. In its traversing stage a single dense group of 
boids (the flock), confined within a non-spanning section of the lattice, slowly flows across the 
system at a consistent non-zero velocity and diffuses. At a later time the reversal mechanism 
begins: this is initiated by a fluctuation near the front of the flock causing a local reversal in 
momentum. The fluctuation moves against the direction of the flock and grows in density and 
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Figure 2: Diagram in the k-/3/(1 + (3) plane showing the regimes present in a system with 
L = 128, (p) — 4, rj — 0.05 with alignment and centring effects only (C — > oo). The dipoles, 
alternating flock and homogeneous flock all occupy characteristic regions. Intermediate states 
at the boundaries of these regions are characterised by mixed behaviours: (a) Flocks behave 
for sustained periods as both homogeneous flocks and alternating flocks, (b) dipoles begin 
to traverse large distances in single motions. The model of ^J] is case (c), a homogeneous 
flock for this parameter set and system size. 



9 



momentum until it passes entirely through the flock inverting its velocity. At this time the 
traversing behaviour of a flock with opposite velocity is established. Initiation of reversals 
can occur either internally or by collision with boids external to the flock. While not all 
large momentum fluctuations initiate a full reversal, they occur often enough to reverse the 
flock before it spreads out to occupy the whole lattice through diffusion. The reversals occur 
stochastically with no fixed period but there is a well defined mean time between reversals. 

In the alternating flock instabilities occur at the front edge of the flock. In it was 
established that the reversal timescale for alternating flocks is 0(ln N) where N is the number 
of boids. The argument is that this is the timescale on which the front edge becomes 
susceptible to relatively large momentum fluctuations, due to its low density. In our model 
there is a new form for G, but logarithmic timescales appear to persist for the alternating 
regime as shown in figure Eb- 

As was shown in figure from many initial conditions a homogeneous flock develops 
after an initial transient alternating flock. However many homogeneous flocks are temporarily 
unstable towards an alternating flock profile even after a long time in the homogeneous 
regime. Such a process occurs for the homogeneous flock of figure after 58000 timesteps. 
In the breakdown process a fluctuation in momentum, initially confined to a few sites, creates 
a disturbance which moves through the flock in a manner similar to the reversal mechanism 
in alternating flocks. Figures Hk,b,c show the progression of the fluctuation through the 
flock to a point where global velocity is completely inverted. Before the flock returns to a 
homogeneous flocking regime there are several further reversals. The localisation of the boids 
and high reversal rate during this period characterise a transient alternating flock state (see 
figures |Ul,e). 

It appears that the temporary failure of the homogeneous flock arises out of a large local 
fluctuation in momentum. During flock reversal the initial fluctuation acquires boids from 
the oncoming flock and continues to grow in momentum and stability. Although it appears 
that the size and type of disturbance required might vary for different systems, the rate 
at which such disturbances occur over the system should be proportional to the number 
of nucleation sites. Within the homogeneous flock, there is a homogeneous density so the 
number of nucleation sites is proportional to L, thus we expect the breakdown time for the 
flock, Tb to behave as 1/L. In figure^ we take the mean time a homogeneous flock of positive 
momentum first reverses and attains a negative global velocity as the flock breakdown time. 
We fitted the data by Tg « j + AL~ B , where we expect B to be approximately 1. The 
constant h is our approximation to the time for the fluctuation to pass through half the flock. 
Such a prediction appears to closely match computer simulations as illustrated in Figure Et- 

The previous paragraph implies that for large systems the homogeneous flock will always 
become unstable due to nucleation of momentum fluctuations. To quantify this, consider 
how the breakdown of the homogeneous flock to an alternating regime is initiated by a 
sequence of successive microscopic flips against the preferred direction sufficient to invert 
the momentum. The number of flips required for this to occur will be proportional to the 
local density, (p), and number of sites over which G is determined (range, r = 3 in our 
simulations). Given that a flip occurs with probability Pp = — 1— -, flock breakdown will be 
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Figure 3: a) Breakdown times for homogeneous flocks and b) reversal times for alternating 
flocks vs InL where L is the system size. Data is averaged over 100 and 1000 breakdowns 
and reversals respectively, (a) HI: [rj = 0.0413,/? = 1, re = 0, (p) = 10], H2: [rj = 0.125, (5 = 
oo, re = 0, (p) = 8]. The straightline is a least squares fit to T = j + AL~ B (see text). The 
exponent is B = 1.03 for H2, and 0.89 for HI. (b) Al: [rj = 0.005, (3 = 2, re = 0, (p) = 1], 
A2: [?7 = 0.02, /3 = 8, re = 0.25, (p) — 1, C = oo]. The reversal times in alternating flocks are 
logarithmic and the straight lines are least squares fits to the data. 
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Figure 4: The Homogeneous flock of \T]p undergoes a rapid change of behaviour. 
(a,b,c) Density and G are plotted as a function of position at fixed times, a) t =58014: 
values are consistent with a homogeneous flock throughout except near site 32 where a 
localised momentum fluctuation exists. (b)t =58046: The interface between the two opposing 
domains created moves leftwards and the fluctuation gathers momentum. (c)t = 58078: 
Within a short space of time the entire flock is reversed, and the density and G profiles are 
characteristic of an alternating flock. Following this event an alternating flock is temporarily 
established, before an homogeneous flock reemyerges. 

(d,e) A measure of the spatial distribution of boids and the global velocity plotted against 
time from before the breakdown, through the temporary alternating flock until the emergence 
of a homogeneous flock. 



initiated with probability ~ LP F r in any timestep. Therefore one expects the typical time 
r for the flock to flip to be r ~ L~ 1 Pp kr ^ p '. This implies that for the homogeneous flock to be 
stable over times r ~ L a where a > 0, a density (p) ~ In L would be required. We conclude 
that for finite range and density no homogeneous flocks will be stable in the thermodynamic 
limit (L — > oo). However, as we have seen in Figure El on intermediate size systems there are 
well defined regions of parameter space where homogeneous flocks do exist. 

Finally we examine the impact of varying the strength of centring and separation on 
flocking states. The alternating flock persists even when centring is a relatively strong effect, 
whereas the homogeneous flock typically becomes unstable towards an alternating flock under 
such conditions. Thus in the phase diagram in figure El for large enough j3, increasing k leads 
to transitions from the homogeneous to the alternating flock and thence to the dipole state. 
The phase diagram with the inclusion of separation (finite C) is qualitatively similar. Perhaps 
surprisingly, decreasing C does not substantially affect the alternating flock. 



3.3 Dipole Regimes 
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Figure 5: [rj — 0.1, (3 — 4, k — 0.75, C — > oo). Density Profiles (upper panels) and profiles of 
the expected updated velocity, G, (lower panels) after 50000 timesteps. (a) At high density 
(p) = 2000 domains, with exponentially decaying densities and fixed velocity, bounded by 
discontinuities in G are clear, (b) At density (p) = 1 dipoles occupy a small portion of the 
lattice, but contain the majority of boids. There are large portions devoid of boids with 
G = 0. 



The dipole regime is characterised by single or multiple localised structures. These may 
be stationary or slowly moving. The density is maximal at the centre of the dipole and the 
density profile is symmetric about the centre. However, to the left of the centre G is positive 
but to the right of the centre G is negative. Thus the structure traps particles in such a way 
that the net flux of boids across the centre of the dipole is zero. 
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In figure|3]we illustrate two systems, one at high density and the other at densities typical 
of our other simulations (approximately 1). In the case of high density it can be seen that the 
system is dominated by adjacent large dipoles. Each of these dipoles has an exponentially 
decaying density to a boundary with another dipole, and a steady value for G trapping boids. 
In this system there is slow redistribution of boids across boundaries towards larger dipoles. 
In the second system density is lower and the exponential density profiles are not obvious; 
instead there is one large dipole and several smaller dipoles. The dipoles span a very small 
fraction of space and are not in direct contact with others. The gaps between the dipoles 
are essentially homogeneous regions of very low density but small dipoles (of size 2-5 boids 
in the figure) regularly form and evaporate away. Occasionally a small dipole arising out of 
the low density region can grow to supersede larger dipoles. 

In figure IH1 details of a typical coarsening process, whereby small dipoles are eliminated 
and large dipoles grow, is illustrated. Figure El shows the development of a system, after 
an initial transient period when many dipoles form, through an intricate coarsening process. 
Figure shows the system coarsening into two large dipoles and eventually (at a later 
time not included in [H^.) one large dipole emerges. This may be seen from figures ^p,c 
which plot the statistics of the sizes of the dipoles remaining in the system. The regions in 
between the dipoles is a low density domain containing independent boids which combine 
sporadically to form small transient dipoles. In the fully coarsened system there is a relatively 
small probability of a boid being elsewhere than the single dominant dipole — all dipoles of 
intermediate size disappear. 

In the low density regime (as in|Hb) the exchange of boids between dipoles occurs through 
the homogeneous low density domains. Dipoles may lose boids through stochastic fluctu- 
ations at their boundaries. If two adjacent dipoles have such noise-induced boid loss at 
different rates, there will be a net flux across this domain towards the more stable dipole. 
Numerically we have observed that in fact smaller dipoles have a higher loss rate than larger 
ones. Therefore the rate of boid loss depends on the fine details of the dipole profile which 
depend on the dipole size. (An assumption of a pure exponential density profile for all dipoles 
would imply no dependence of the rate of boid loss on dipole size.) Other effects might also 
be important, for example the greater mobility of smaller dipoles. 

Dipoles are created and sustained where centring is the dominant effect. Within dipole 
systems increasing alignment (decreasing k) is most noticable in the resulting greater mobility 
or wandering\2'i] of dipoles (and independent boids); a significant process since it hastens 
the coarsening of the system, especially in the initial stages. Wandering can occur in any 
size dipole but is most prominent in smaller dipoles. Typically wandering is initiated by a 
net flux of boids across the centre of the dipole causing the centre to shift. This motion 
can be induced or enhanced by alignment, since the matching of velocities in such a process 
favours the collective motion. With alignment sufficiently strong, this process can be self- 
reinforcing to the extent that a sequence of movements or sustained translations (flocking) 
become possible. 

Dipole systems are characterised by their central region of high density. However, we 
might expect such a profile to be suppressed by separation. The inclusion of separation in 
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Figure 6: System with same parameter set as figure on a larger system of L = 1024. (a) 
A typical space-time plot starting from random initial conditions. After an initial transient 
period when the dipoles form, an intricate coarsening process takes place wherein small 
dipole eventually dissolve and large dipoles grow. (b,c) Statistical data on the remaining 
dipoles as a function of time averaged from 100 runs similar to (a), (b) The number of dipoles 
remaining of size greater than the values indicated in the legend: Eventually the number 
of dipoles of intermediate size decreases to zero (curves B,C,D,E converge to 1); After 10 7 
timesteps a single large dipole emerges alongside approximately 2 small transient dipoles. 
(c) The probability that a boid is contained within a dipole greater than the sizes indicated 
in the legend as a function of time: Only a small number of boids ever exist outside dipoles 
(A) and after sufficient time the ma ioritv of boids are contained in a, single larere dinole CE). 
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Figure 7: [rj = 0.02, /3 = 8, k = 0.85, (p) — 1,C — 20] This system contains two large 
fixed-density dipoles alongside small dipoles and solitary boids. In small dipoles centring is 
the dominant effect and wandering is limited, however in the two fixed-density dipoles large 
translations are observed. See text for explanation of the term fixed-density dipole. 

dipole systems leads to a dramatic alteration of behaviour: where alignment is a small effect, 
and (p) < C < oo, a system of fixed- density dipoles can be observed. Fixed-density dipoles 
have reduced density at the centre, but the density is constant within an extended central 
region and only near the boundaries does one see the sharply decaying density profile, thus 
the density profile has a flat top. Within fixed-density dipoles alignment plays a much more 
important role than in comparably sized dipoles in the absence of separation, since at the 
centre of fixed-density dipoles the separation and centring tendencies tend to cancel allowing 
alignment to dominate. The relative stability of large and small fixed-density dipoles is left 
as an open question. Fixed-density dipoles can, for example, split from their centre and 
wander more substantially. These processes are appreciable in figure [7| where the mobility 
of the large fixed-density dipoles is second only to that of solitary boids. 

4 Mean Field Theory 

In order to understand some of the regimes observed numerically above, we now develop, 
from a mean field treatment of the full dynamics, a continuum model in which the density p 
and neighbourhood velocity v are continuous functions. In so doing we are able to identify a 
linearly stable flocking solution comparable to the homogeneous flock and a piecewise con- 
tinuous solution comparable to the dipole regime. Furthermore we show that these solutions 
exist in numerical iterations of the mean-field equations. 
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4.1 Equations 



The equations satisfied by the mean-field density, p, and mean-field neighbourhood velocity, 
v, are derived in appendix A and read 



dp 
dt 
dv 
di 



d . , Id 2 



where 
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dx 
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tanh 3 
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uj 



3p dx 



(15) 
(16) 

(17) 



Equation ()15|) is rather easy to understand: note that, since G is the mean velocity of 
boids, pG is a mass current, thus the equation is a continuity equation with a current and 
a diffusive term. The equation for the neighbourhood velocity v f!16|) . on the other hand is 
more complicated. 



4.2 Steady State Solutions and Linear Stability 

We now look for steady state solutions of the above equations, for which we set the lhs of 
(I15I16|) to zero. Further, we look for homogeneous solutions where spatial derivatives are 
zero. The global density may be freely chosen, in these solutions p(x,t) = (p), whereas 
velocities are determined by (fTB^) as v = G(v, 0) which yields 

V = {1 ~ 27]) ^ • (19) 

Thus, there are either one or three solutions to the mean field theory which are homogeneous 
in density. 

For all parameter sets one solution to this is v = 0, the disordered solution, whilst 
for certain values of (3, k and rj two flocking solutions v = v± exist, corresponding to the 
homogeneous flock. The existence of the flocking solutions is determined by the gradient in 
G with respect to v at the origin. Since G is a concave function of \v\ a flocking solution of 
(HHJ) requires that G'(0, 0) > 1. Thus the existence of this flocking solution is dependent on 
the following criterion being satisfied: 

tanh 8 , . 

< (1-277)(1-k). (20) 

This is possible in systems with low rj, but strong alignment (so that 8(1 — k) is large). 
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4.3 Piecewise Dipole Solution 



We show in this section how it is possible to determine a further steady-state solution of 
(jl5H6|) where the density is an exponentially increasing or decreasing function of x and v 
is a positive or negative constant respectively. Domains of these solutions may be pieced 
together so that the density is continuous to form dipole structures. 

First, we assume a fixed value of u implying by definition (J 18)) a density profile 



p(x,t) 



Aexp 



3cux 



(21) 



where A is a constant. Next we propose a solution whereby within any domain G and hence 
velocity are fixed. Then, from (|THj) we have 







3 9 
-G(w,v)-up(x) + -uj 2 p(x) 



which implies 



G(u,v) 



3u 

T 



(22) 
(23) 



The velocity v must similarly be steady and ()16j) implies 

: 



3 3 9 9 

—v + G(lu, v) H — u)G(u, v)v to uj 2 v H — u 2 G(u, v) 

2 2 8 8 



which after substitution of G from 



3u 

T 



simplifies to 

2 



2't 



(24) 



The condition required for a non-homogeneous solution is the consistency of f!23l) and 
(I17|) . eliminating G we obtain 



4(1 -2r?) , , , 

lu = — tanh { (3 

3 tanh p 



3 



3oo 

T + 2(l 



3oj V 



(25) 



Alongside the homogeneous density solutions, = 0, there are potentially dipole solutions 
of (|25j) w± 7^ 0. Correspondingly u takes values u (cj ± ) via ()24p . One of these has rightward 
G and exponentially increasing density, the other a leftward G and exponentially decreasing 
density. A dipole is a localised structure consisting of a left hand region of exponentially 
increasing density and a right hand region of decreasing density. To satisfy periodic boundary 
conditions it is necessary that the system consists of complementary domains. Consider a 
boundary between the dipole domains: the two exponentially decaying or exponentially 
growing domains meet with equal density at such a point. The density gradient and other 
qualities are locally antisymmetric about the boundary; hence the net flux of boids will be 
zero and the boundary stationary (this would not be true for a dipole domain in contact with 
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a homogeneous domain). Within each domain the presence of a boundary is not felt, and 
so domains of differing width and height can neighbour one another. Hence non-symmetric 
dipoles of differing sizes can be neighbours. Our solutions are therefore the set of coupled 
non-overlapping dipole systems (domain pairs i), each defined uniquely by the set of dipole 
centre positions, heights and widths. 

These solutions are comparable to the full stochastic system at very high density (see 
figure EK)- Note that the larger dipoles in the sparser system (figure Eb) also approximate a 
pair of domains with the same fixed G and uj. 

In numerical iterations of the mean field equations (see next subsection) the dipole struc- 
tures were rounded at the cusps which form at domain boundaries. One expects this rounding 
and the physics at the cusps to be described by higher order terms in the lattice spacing a 
which have been ignored (see Appendix A). 

4.4 Numerical Analysis 

Numerical iteration of the dynamical mean field equations (jl5H6|) allows further analysis 
of these states and comparison to full system dynamics. The iteration was carried out by 
rediscretising (|15|l(jj) onto a lattice. We found the flocking (uj = 0,v±), non-flocking {uj = 
0, v = 0) and dipole solutions (uj±,v(uj±)) to be present, in addition to systems comparable 
to the fixed-density dipole. The alternating flock however was never observed, nor were 
systems comparable to sparser dipole regimes. 

Where a flocking solution is predicted, and not made unstable by the presence of centring 
(see section I4.6J1 a homogeneous flock with the anticipated global velocity (jTUjl was always 
observed to emerge from any initial conditions. 

In dipole regimes we found that the system iterated to a dipole system consisting of 
many dipoles. However the mean field dipole system does not exhibit a further coarsening 
process comparable to the full stochastic model. This is likely due to the absence of noise 
driven fluctuations which drive the coarsening process. On the other hand, where the initial 
boundary between domains does not fall perfectly on a lattice space boundary or centre a 
very slow redistribution of boids does occur in the mean-field system, but apparently at only 
an exponentially slow rate. This redistribution was not observed to cause the disappearance 
of any dipoles, nor did it appear to favour transfer to larger dipoles (is the case in the full 
stochastic system). 

4.5 Linear Stability of Mean Field Solutions 

We have developed a linear stability analysis of the disordered and flocking solutions to 
the mean field equations. This is done by linearising (J 151 computing the evolution of 
the Fourier components then checking whether positive (unstable) eigenvalues exist. As the 
details are rather involved we do not present them here but quote some results. In the case 
k = 0, when the flocking solutions exist they are linearly stable and the disordered solution is 
unstable; when the flocking solution does not exist the disordered solution is stable. However, 
for sufficiently large k an instability will occur in both flocking and disordered solutions. 
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This instability is expected to be with respect to stable dipole solutions, although we did 
not explicitly check the stability of dipole solutions. 

Thus the equality in (|2U|) defines a mean-field phase boundary between the disordered 
and possible flocking states. Comparing to Figure 121 this curve is in approximate agreement 
to the boundary between the disordered and alternating flock regimes. However for other 
parameter sets the boundary is less accurately predicted. Generally the boundary becomes 
increasingly more accurate as the density increases. 

As well as the disordered/flocking phase boundary, the instability of the flocking solutions 
at large k is also consistent with Figure El where for large k dipole solutions appear. The 
only discrepancy remains the absence of alternating flocking solutions within the mean field 
theory: in regions where one sees an alternating flock in the full stochastic system the mean 
field theory has a homogeneous flock solution. 

4.6 Absence of Alternating Flock and Stability of Mean Field So- 
lutions 

As noted above the alternating flock was not expressed within the mean field theory, even 
in the time dependent form (fH)l fTo^ . As we shall now discuss, the mean field theory is 
not capable of describing the alternating flock which is driven by fluctuations. However 
the flocking solution of the mean field theory does appear to correctly predict the velocity 
of the alternating flock between reversals i.e. the solutions of (fTT)|) proves to be in good 
agreement with transient results from the full stochastic system. We also note that we found 
numerically that the mean field theory correctly describes the evolution of the profile of an 
alternating flock as it spreads between reversals. However, whereas an alternating flock in 
the full stochastic system will inevitably reverse due to some fluctuation in momentum (as 
described section 13 .2|) the mean-field profile will keep evolving until the homogeneous flock 
is attained. Similarly, due to the absence of fluctuations the homogeneous flock will not 
reverse for any system size. 

Clearly, to reproduce the alternating flock, a noise term should be added to the mean 
field equations (fTol [TBJ) thus turning them into Langevin equations. This would be done, 
following the proposal of [E1HEI, by adding a noise term to the velocity equation (|TB|). The 
equations (fTol [TB|) would then be similar to those of [121 HE] although our velocity equation, 
which we derived from the microscopic dynamics, is rather more complicated and contains 
additional terms. 

5 Discussion 

In this work we have extended the flocking model of ^H] to include all three of Reynolds' three 
effects — alignment, centring and separation — implemented in combination or independently. 
Within this model we demonstrated the robustness of the alternating flock highlighted in 
[TU] to the addition of other Reynolds' effects. Also we showed the existence of two new 
regimes: the homogeneous flock and dipole structures. These are consistent with results 
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from a mean field treatment which, in particular, correctly predicts the form of the dipole 
solutions. Furthermore, we investigated the coarsening process in the low density dipole 
regime. 

We also derived from microscopic considerations the proposed form for the function G ©• 
The algorithm presented in section 12.31 is not only intuitive but significantly is a fixed time 
algorithm, by which it is meant that the algorithm would involve sampling a fixed number 
of neighbouring boids whatever the density. Such algorithms are aspired to in flocking 
applications [3J E] and are presumed to underly natural boid decision making 7 . Several 
ways to extend the algorithm could be to include a density, velocity or historical dependence 
in determining the sample, anisotropic spatial sampling or considering a direction selection 
rule other than majority. For example boids could require unanimity[20 amongst their 
observed neighbours to believe that any orientation is a worthy choice. 

Reynolds' three effects were defined as required behaviours for simulating real flocks. In 
this paper we have shown several states which may characterise real systems. Firstly the 
homogeneous flock, which has many analogies in two or three dimensions, can be realised in 
some systems approximating one dimension: for example skaters confined to the outside of 
a rink, or fish in a ring shaped tank (as in the Boston Aquarium). The fixed density dipole 
is characterised by sharp edges and fixed internal density, and the ability to wander — either 
slowly, or with the inclusion of alignment, quickly and with a sustained orientation. These 
qualities, as well as as the ability to split from within, could describe a number of biological 
systems. Finally consider the reversal mechanism present in the alternating flock; such sharp 
reversals of direction are characteristic of the manoeuvres seen in many flocks, though they 
are poorly described by many models. 

The numerical results obtained for the full stochastic system have some noteworthy simi- 
larities with previous studies. A flocking state is observed in the one dimensional continuous 
space model of JB] which appears to resemble an alternating flock. Dipole-like regimes have 
been observed in [21], where dense states are supported by a centring interaction. As with 
our dipoles the density at the centre of these systems grows rapidly with boid number. This 
effect is suppressed by a hard core repulsion, producing states similar to our fixed density 
dipoles. In [21] two dimensional systems also display wandering and oscillating (circling) 
behaviour. In a two dimensional cellular automaton demonstrates a density-dependent 
state with symmetry-broken velocity. This density dependence arises from a hard capacity 
being placed on the number of boids which can occupy a single site; in effect the flocking 
state is destroyed by separation at high density. 

Our mean field treatment has similarities with the continuum model in JH]- Not only 
are the same states demonstrated, but the terms within the equations are comparable also. 
Special note is made in JH] of a term of the form ^J^Gj^p which is present in our equation 
(I16j) . This term relates to domain competition, and aids in an understanding of the reversal 
mechanism and breakdown of the homogeneous flock. As (jlfij) describes very well the shape 
of a spreading flock, and it is believed that this shape allows certain fluctuations to cause the 
flock reversal mechanism JH] > it might be hoped that the addition of noise into ()15I16|) would 
allow the alternating flock to be further analysed. The coarsening process in low and high 
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density dipole regimes may also be observed by allowing fluctuations to enter the dynamics 
at the boundaries of the dipoles. 

Finally we mention one further variation of the model of ^H] that we have studied |20| . 
To imitate the effects of inertia we modified step (iii) of our dynamics such that boids which 
decide to change direction must proceed through a temporary zero velocity state. Numerical 
simulations and a mean field treatment (which we do not present here) indicate that under 
such a scheme flocking solutions become more prominent and have higher global velocities. 
Although the flipping mechanism proceeds in much the same manner as the two velocity 
case, the rate of reversal for both the alternating and homogeneous flock is reduced This 
is not surprising since more unfavoured flips are required to create momentum fluctuation 
leading to a sustainable reversal of the velocity. 
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A Derivation of Mean Field Equations 

In our mean field theory we shall approximate the density p(x, t) (the number of particles at 
site x and time t) and the momentum <j)(x, t) (the number of right going minus the number 
of left moving particles at site x and time t) by continuous functions and derive the evolution 
of these funcions by ignoring certain correlations. We begin with the dynamics of site x in 
a single update, as defined by p{x,t) and <p(x,t). Consider p(x,t), in a single update this 
can increase if a boid selected from a neighbouring site moves into site x, or can decrease if 
a boid is selected from site x, with an associated probability: 

p(x,t + 6t)= (26) 

p(x, + 1, Probability = ^^W+{x - a,t) + ^p^W. (x + a, t) 
p(x,t) — 1, Probability = 
p(x,t), otherwise. 

We then average over the events occurring between time t and t + 5t 

p(x,t + 5t) = Yl p(x,t + St)P KXtt+St) (27) 

p(x,t+5t) 

and make the mean field approximation of factorising all averages. 

Further we assume p(x, t) and (f>(x, t) vary smoothly on the time and length scales effective 
in a single update. With these assumptions we expand about site x to second order (in lattice 
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spacing a) and about t to first order (-^). 

dp r dp a 2 d 2 p^ rTJT dW + a 2 d 2 W + ^ 

dp a 2 d 2 p, rTJT dW- a 2 d 2 W_, 

+ {p + a lT x + ^ 2}{w ~ + a l^ + Yl^ } 

- P- (28) 

By introducing G © the following equation can be determined. In this equation there 
is a diffusion term arising from the choice of random stochastic updates, and a current term 
controlled by the mass current pG. 

In a similar way one can develop an equation for <fi(x, t) beginning with the full dynamics at 

x, t, 

( 4>{x, t) + 1, Probability = B^^W+{x - a,t) + 
<f){x, t + St) = I t) - 1, Probability = ^p^W^x + a,t) + 
y (j)(x,t), otherwise. 

Taking the expectation value and expanding to second order in a and first order in t: 

d(j> r dp a 2 d 2 p^ rTjr dW+ a 2 d 2 W + ^ 
dt = {f> - a d-x + ^ }[W+ ~ a ^x- + ^^x 2 - } 



dp a 2 d 2 p, ( „ T dW_ a 2 d 2 W_, 



(30) 



Using (JHJ this can be simplified to its final form 



Defining the (neighbourhood) velocity as 

yields using (|29I31|1 a velocity equation: 

dv ad d a 2 d 2 d 2 

We now write G(x,t), the expected update velocity, as a function of v(x, t) and p(x,t) 
only. To simplify we keep only leading order terms in a and obtain 

G=(1 _ 2ty) tanh(/?[(l-^ + H) (33) 
24 



where 



2a dp(x,t) 
3p dx 

Finally, setting a to one in J2SEEJE1 yields (fToUTS|). 
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